Temporally auto-correlated predator attacks structure ecological communities

For species primarily regulated by a common predator, the P* rule of Holt & Lawton (Holt & Lawton, 1993. Am. Nat. 142, 623–645. (doi:10.1086/285561)) predicts that the prey species that supports the highest mean predator density (P*) excludes the other prey species. This prediction is re-examined in the presence of temporal fluctuations in the vital rates of the interacting species including predator attack rates. When the fluctuations in predator attack rates are temporally uncorrelated, the P* rule still holds even when the other vital rates are temporally auto-correlated. However, when temporal auto-correlations in attack rates are positive but not too strong, the prey species can coexist due to the emergence of a positive covariance between predator density and prey vulnerability. This coexistence mechanism is similar to the storage effect for species regulated by a common resource. Negative or strongly positive auto-correlations in attack rates generate a negative covariance between predator density and prey vulnerability and a stochastic priority effect can emerge: with non-zero probability either prey species is excluded. These results highlight how temporally auto-correlated species’ interaction rates impact the structure and dynamics of ecological communities.


Introduction
Predation or resource limitation can regulate populations. When multiple species are regulated by the same limiting factor, long-term coexistence is not expected under equilibrium conditions. Regulation due to a common, limiting resource can result in the R* rule: the species suppressing the resource to the lower equilibrium level excludes other competitors [1][2][3]. Regulation due to a common predator can result in the P* rule: the prey species supporting the higher equilibrium predator density excludes the other prey species [4][5][6]. Yet, many coexisting species share a common resource or a common predator. Understanding mechanisms permitting this coexistence is central to community ecology. One of these coexistence mechanisms, the storage effect for competing species, relies on temporal fluctuations in environmental conditions [7][8][9]. Whether an analogous, fluctuation-dependent mechanism exists for species sharing a common predator is studied here.
Similar to species competing for a common resource, species sharing a common predator can exhibit mutually antagonistic interactions [10]: increasing the density of one prey species leads to an increase in predator density and a resulting increase in predation pressure on the other prey species. Thus, to the uninformed observer, the prey appear to be competing. Empirical support for apparent competition is extensive [4,11,12] and has significant implications for conservation biology [13]. When the shared predator is the primary regulating factor, Holt & Lawton [5] demonstrated that one prey excludes the other via the P* rule. Yet in nature, coexisting species often share common predators. Holt & Lawton [5] found that spatial refuges, resource limitation and donor-controlled predation could help mediate coexistence. However, environmentally driven fluctuations in demographic rates did not promote coexistence [5,14]. These studies, however, assumed environmental fluctuations are temporally uncorrelated.
By contrast, environmental fluctuations are known, both theoretically and empirically, to mediate coexistence for species competing for a common resource. In a series of influential papers [7,8,15,16], Chesson identified two fluctuationdependent coexistence mechanisms: nonlinear averaging and the storage effect. Empirical support for these mechanisms exist in a diversity of systems [9,[17][18][19][20][21][22]. A key ingredient for the storage effect is a positive covariance between favourable environmental conditions and species' densities. Temporal auto-correlations, which are commonly observed in environmental factors [23,24], can generate this positive covariance [25].
Here, temporally auto-correlated fluctuations in demographic rates are shown to mediate coexistence of prey species primarily regulated by a predator and to generate stochastic priority effects. To derive these conclusions, stochastic models of predator-prey interactions are studied using a mixture of analytic and numerical methods.

Model and methods
Following Nicholson & Bailey [26] and Holt & Lawton [5], the model considers two prey species with densities N 1 , N 2 that are regulated by a common predator with density P. In the absence of predation, the density of prey i increases by a factor, its finite rate of increase R i (t), in generation t. Individuals of prey i escape predation with probability exp (−a i P) where a i is the attack rate on prey i. Captured individuals of prey i are converted to c i predators. To ensure population regulation, predators immigrate at rate I > 0 [5,14]. Allowing for fluctuations in the demographic rates, the model becomes c i ðtÞN i ðtÞð1 À expðÀa i ðtÞPðtÞÞÞ þ IðtÞ: ð2:1Þ Consistent with meteorological models [27,28], fluctuations in logarithmic demographic rates are modelled as a firstorder auto-regressive processes (see electronic supplementary material, (A1) in appendix). For example, the log attack rates ln a i (t) are characterized by their means ln a i , their variances s 2 i ¼ Var½ln a i ðtÞ, their temporal auto-correlation ρ = Cor[ln a i (t), ln a i (t + 1)] and their cross-correlation τ = Cor[ln a 1 (t), ln a 2 (t)]. For the numerical simulations, the auto-regressive processes are Gaussian, i.e. the attack rates are lognormally distributed.
The dynamics of (2.1) are explored using analytical and numerical methods. The analytical methods rely on the invasion growth rates (IGRs) of the prey that correspond to the average growth rates when the prey species becomes rare [7,[29][30][31][32]. When both IGRs for the prey are positive, both species increase when rare and coexist. When both IGRs are negative, there is a stochastic priority effect, i.e. with non-zero probability either prey species is excluded.
When the IGRs have opposite signs, the species with the positive IGR may exclude the other species. Analytical approximations for these IGRs are derived for small environmental fluctuations and computed numerically using R (details in electronic supplementary material, appendix).

Results
The invasion growth rates (IGRs) of the prey species are defined by assuming one species, say species j, is common (the resident) and the other is infinitesimally rare (the invader), say species i ≠ j. The resident prey and predator species are assumed to coexist (see electronic supplementary material, condition (A2) in appendix) and have reached a stationary distribution. Let P j (t) be the predator densities at this stationary state. At stationarity, the average intrinsic growth rate of prey j equals the average predation rate: ð3:1Þ When introduced at (infinitesimally) low densities, the IGR of prey i equals the difference between its average intrinsic growth rate and its average predation rate (see electronic supplementary material, (A5) in appendix) We use equations (3.1) and (3.2) to show the P* rule holds when a i (t) are temporally uncorrelated even if the other demographic rates are auto-correlated. Then we show how auto-correlations in a i (t) generate alternative ecological outcomes.
(a) Temporally uncorrelated attacks and the P*-rule If the attack rates are temporally uncorrelated (ρ = 0), then the average attack rate of the resident prey equals the product of the average attack rate and the average predator density: a j P j ¼ a j P j . Consequently, equations (3.1) and (3.2) imply prey i's IGR is proportional to the difference in the average predator densities supported by prey j and prey i, respectively: Hence, if prey 1 supports the higher predator density (P 1 . P 2 ), then r 1 > 0 and r 2 < 0 and prey 1 excludes prey 2. The opposite conclusion holds if prey 2 supports the higher predator density (P 2 . P 1 ).

(b) The resident and invader attack covariances
Unlike the temporally uncorrelated environments, autocorrelated predator attack rates generate a covariance Cov[a j , P j ] between the predator attack rates a j (t) and the predator densities P j (t) when prey j is the resident. This resident attack covariance depends on the auto-correlation ρ of a j (t) in a nonlinear fashion ( figure 1a). An analytical approximation (see electronic supplementary material, equation (A8) in appendix) shows the resident attack covariance is negative when either the auto-correlation ρ is negative (ρ < 0) or ρ is greater than the reciprocal of the prey's mean finite rate of increase (1=R j , r , 1). Only for positive, but not too positive auto-correlations (0 , r , 1=R j ) is the resident attack covariance positive. Figure 1b,c provides a graphical representation of this nonlinear relationship. When the predator's attack rate increases, there is a short-term increase in the predator's density. However, the continual decrease in the resident prey's density ultimately results in a reduction in the predator density. When the auto-correlation is sufficiently positive, increases or decreases in attack rates persist for a long time and, thereby, generate a negative resident attack covariance. Less positive auto-correlations or negative auto-correlations play out on shorter timescales and generate positive or negative covariances, respectively.
Temporally auto-correlated attack rates on prey i (the invader) also generate a covariance between their attack rates a i (t) and the predator densities P j (t) (figure 1d). The sign of this invader attack covariance Cov[a i , P j ] is determined by the cross-correlation τ = Corr[ln a i , ln a j ] between the attack rates. When this cross-correlation is positive, the invader and resident attack covariances have the same sign (two darkest lines in figure 1d). When this cross-correlation is negative, these two covariances have opposite signs (two lightest lines in figure 1d).

(c) Auto-correlated attack rates alter ecological outcomes
Owing to their effect on the attack covariances, auto-correlated attack rates can alter ecological outcomes ( figure 2a,b). This impact is best understood when the prey only differ in the timing of predator attacks (i.e. R 1 = R 2 , ln a 1 ¼ ln a 2 , s 2 1 ¼ s 2 2 but τ = Corr[ln a i , ln a j ] < 1). Then the IGR of prey i equals the difference between the resident and invader attack covariances: Whenever the prey species experience differential predation in time (τ < 1), the sign of the IGR r i is determined by the resident attack covariance (figures 1a,d and 2b). Hence, if the autocorrelation ρ in the attack rates is positive but not too positive (0 , r , 1=R, dashed lines in figure 2c), the IGRs are positive for both species and the prey coexist ( figure 2a). Alternatively, if the auto-correlation ρ is negative or too positive (ρ < 0 or r . 1=R), then IGRs are negative for both species and the prey exhibit a stochastic priority effect (figure 2b). More generally, when the prey species differ in their intrinsic fitness R i and differ in their mean attack rates, the invasion growth rate of the prey i depends on these differences: The first term log R i À log R j corresponds to the difference in the mean intrinsic rates of growth of the two prey species. When prey species i has a larger mean intrinsic rate of growth, this term is positive, otherwise it is negative. The second term, ða j À a i Þ P j , is proportional to the difference in the mean attack rates. When the common prey species is less vulnerable on average (i.e. a i , a j ), this term is positive, otherwise it is negative. The final pair of terms, the difference in the resident and invader attack covariances, is equivalent to (3.4). Hence, differences in the prey's mean intrinsic rates of growth or mean vulnerability to predation can alter the effects of temporally auto-correlated attack rates. For example, large differences in mean attack rates can lead to exclusion despite the attack covariances helping increase IGRs.

Discussion
For species primarily regulated by a common predator, coexistence is not expected under equilibrium conditions: the prey supporting the higher predator density can exclude the other via apparent competition [5,14]. Similar to the storage effect for competing species [7,33], we found that environmental fluctuations impacting prey specific attack rates can modify this ecological outcome. Two conditions are necessary for these alternative outcomes. First, fluctuating environmental conditions must differentially impact the predator attack rates on the different prey species. This condition is equivalent to 'species-specific responses to environmental conditions' required for the storage effect [7]. The second condition requires a non-zero, within-generation covariance between predator attack rates and predator density. This Dashed lines correspond to analytical predictions of where this covariance vanishes. This nonlinearity stems from the short-term versus long-term effects of an increase in the predator attack rate on the predator (b) and prey (c) densities. In (d ), the invader attack covariance Cov[a i , P j ] is plotted for different cross-correlations τ. Parameters: attack covariance is analogous to the 'environment-competition covariance' of the storage effect for competing species [7]. When positive, the attack covariance results in relatively lower predation rates on prey that become rare and, thereby, facilitates their recovery from low densities. Hence, coexistence is more likely. When the attack covariance is negative, predation rates are relatively higher on prey that become rare resulting in a stochastic priority effect.
There is empirical evidence that suggests both conditions are likely to occur in nature. For the first condition, differences in prey vulnerability to predation provide multiple pathways for generating asynchronous attack rates among multiple prey [34]. These pathways include differences in micro-habitat and refuge availability [35,36], environmental stressors [37], phenology [38] and morphology and behaviour [39][40][41][42]. For the second condition, temporal auto-correlations are ubiquitous in environmental factors that drive these pathways [23,24,43] and, as shown here, can generate covariances between attack rates and predator densities.
We demonstrate that both the sign and magnitude of the temporal auto-correlations in attack rates determine the sign of the attack covariance. For positive auto-correlations, the sign of this covariance depends on the timescale at which the fluctuations occur. When temporal autocorrelations are weak, fluctuations occur on shorter timescales, generate a positive attack covariance, and promote coexistence. By contrast, when temporal auto-correlations are strong, fluctuations occur over longer timescales, generate a negative attack covariance, and promote stochastic priority effects. Similarly, for continuous-time models of species competing for a common resource, Li & Chesson [44] found that fast resource depletion generates a positive environmentcompetition covariance and, thereby, can promote coexistence. This positive environment-competition covariance arises from consumer attack rates being positively auto-correlated at the timescale of the resource depletion.
The work presented here and earlier work [25,44] highlight that covariances between species densities and per capita species interaction rates can fundamentally alter the composition and dynamics of ecological communities. These covariances can be driven by the sign and magnitude of auto-correlated fluctuations in environmental conditions. Importantly, the magnitude of these auto-correlations can lead to different ecological outcomes due to differences in transient versus long-term responses of species to changing interaction rates [45,46]. Understanding how these effects combine across multiple species, how they interact with other coexistence mechanisms [33], and how they are impacted by demographic stochasticity [47][48][49] provide significant challenges for future work.  Figure 2. Auto-correlated attack rates alter ecological outcomes. In (a), the dynamics of coexisting prey species and (b) two realizations of the dynamics of prey species exhibiting a stochastic priority effect. In (c), the invasion growth rates r i for both prey species when rare as a function of the temporal auto-correlation ρ in attack rates. Different lines correspond to different cross-correlations τ = −1, 0, 1 in the attack rates. Dashed lines correspond to where the analytic approximation of IGRs r i are zero. Parameters: R 1 ¼ R 2 ¼ 2, I ¼ 10, ln a 1 ¼ ln a 2 ¼ ln 0:04, s 2 i ¼ 0:25, c 1 ¼ c 2 ¼ 1 and τ = 0.5 for panels (a,b).